!p.multi=[0,2,1]
!p.charsize=1.5
l=188
rr=8.51904e8         ;m
rmin=0.15
rmax=0.97
rds=rmin*rr     ;m
rdsi=1./rds
GG=6.674e-11   ;N(m/kg)^2
M=0.7*1.988e30     ;kg
dz=(rmax-rmin)*rr/double(l)
r=rds+indgen(l)*dz
g=GG*M/r^2     ;m/s^2
g0=g[0]
z = indgen(l-1)*dz

; bottom values

rho00= 3.29
p00=   1.71889e+13
T00=   2.97382e+06
th00=  1.14*T00       ;2.06008e+06

Rg=13732.
cp=2.5*rg
cap=rg/cp
capi=1./cap
ss=g0/(cp*th00)

rho = rho00*(1.-(ss*z)/(1.+z*rdsi))^(capi-1)
;pm0=p00*(1.-cap*z*sdi)**capi ; CAUTION: max. alt. exists
;tm0=T00*(1.-cap*z*sdi)       ; CAUTION: max. alt. exists

plot,r/rr,rho

st=-3./(th00*(l-1)*dz)
the=th00*(1.+ st*z)
plot,r/rr,the,yr=[min(the),max(the)]

!p.multi=0
end

